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Abstract 

We show that at the onset of a cyclic fold bifurcation, a birhythmic medium composed of gly- 
colytic oscillators displays turbulent dynamics. By computing the largest Lyapunov exponent, the 
spatial correlation function, and the average transient lifetime, we classify it as a weak turbulence 
with transient nature. Virtual heterogeneities generating unstable fast oscillations are the mecha- 
nism of the transient turbulence. In the presence of wavenumber instability, unstable oscillations 
can be reinjected leading to stationary turbulence. We also find similar turbulence in a cell cycle 
model. These findings suggest that weak turbulence may be universal in biochemical birhythmic 
media exhibiting cyclic fold bifurcations. 
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I. INTRODUCTION 



In studies of chemical turbulence in reaction diffusion systems near a Hopf bifurcation, 
a reduction of the model to the complex Ginzburg Landau equation(CGLE) is very useful 



First, it allows to define a parameter set in the model leading to turbulence without 

n 



DD 

carrying out simulations |3[. Secondly, the detailed knowledge of the CGLE's dynamics can 
be very helpful □ □□□ , because mathematical models from different disciplines displaying 
dynamics near a Hopf bifurcation obey the same qualitative dynamics of the CGLE [8[ . 

However, the CGLE alone is insufficient for a qualitative description of realistic models 
in a neighborhood of a Hopf bifurcation, when other bifurcations occur nearby 
For example, near a supercritical Hopf bifurcation point, another stable limit cycle may 
exist, so that, depending on initial conditions, oscillations with two different frequencies and 
amplitudes are possible. Such a situation called birhythmicity is a characteristic feature of 
a number of well known models in biochemical oscillations [la . Il^ . For these systems, the 



CGLE cannot be used without appropriate modifications. Often, the best way to approach 
these problems is simulations of the original models jiol 

To the best of our knowledge, little is known about turbulence in birhythmic media. 
Intuitively, in a regime of a strong wavenumber instability, birhythmicity should not be a 
factor. Therefore, turbulence in homogeneous birhythmic media and in coupled limit cycle 
oscillators should have similar characteristics. In the absence of wavenumber instability, 
high frequency oscillations are supposed to suppress slow oscillations and restore uniform 
oscillations. But at the onset of a cyclic fold (CF) bifurcation in birhythmic media of a 
biochemical origin, high frequency oscillations may be unstable. Thus, a complete suppres- 
sion of slow oscillations may not be achieved in these systems. On the contrary, if unstable 
oscillations emerge persistently, complex spatiotemporal motions are possible. 

The goal of this work is to show that near cyclic fold bifurcations in birhythmic media, 
virtual heterogeneities creating unstable oscillations can lead to a peculiar turbulence, inter- 
mittency of small and large amplitude oscillations. We will first compute complex spatiotem- 
poral behavior in a birhythmic medium composed of glycolytic oscillators. By calculating 
the maximal Lyapunov exponent, the spatial correlation function, and the average tran- 
sient lifetime, we will provide evidence that this behavior is weak transient turbulence. In 
the presence of wavenumber instability transient turbulence may become stationary. Math- 
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ematically, the instability of the faster oscillations is a result of a CF bifurcation driven 
by the terms representing enzymatic regulations, suggesting that weak turbulence may be 
common in biochemical birhythmic media exhibiting CF bifurcations. As further evidence, 
we demonstrate weak turbulence in a cell cycle model. A biological system where weak 
turbulence might possibly be found is presented in the closing section. 



II. A BIRHYTHMIC MEDIUM OF GLYCOLYTIC OSCILLATORS 

Let us introduce a birythmic medium composed of glycolytic oscillators, 



da Oi'i n „ . 

— = Qa<P- k sl - Rn + T + £> 7 A 7) (2) 

q(l + a)(l + 7 ) 2 
9 L + (l + a) 2 (l + 7 ) 2 ' 

In Eqn. (1-2), a and 7 represent dimensionless substrate and product concentrations 
of glycolytic reactions, K, n, u, <Ji, a, k s , L and Q are parameters. For convenience, we 
assume Q = 1 throughout this paper. D a , D 1 are diffusion constants for the substrate 
and product. Our units of time and space are sec and cm, respectfully. When D a = 0, 
D 7 = and er^ = 0, Eqn. (1-2) is called the glycolytic oscillator. The term y ^^ n represents 
substrate recycling that drives birhythmicity. Recently, in Ref. Eqn. (1-2) were shown 
to support multiple wave fronts. Our concern in this paper is a different parameter region 
where irregular spatiotemporal motions develop. 

A phase plane analysis of Eqn. (1-2) shows that the mechanism of birhythmicity is 
two regions of negative slope in the product nullcline Q]. A convenient way to illustrate 
birhythmicity is a bifurcation diagram. We used a well known software package, AUTO [16], 
for bifurcation analysis of the local mode^-Do, = -D 7 = in Eqn. (1-2)). Solid lines in Fig. 
1 show stable steady states, dashed lines show unstable steady states. Stable limit cycles 
are shown by filled symbols, unstable limit cycles by open circles. Filled circles represent 
small amplitude oscillations with high frequencies. Large amplitude oscillations with low 
frequencies are shown by filled diamonds in Fig. 1. A Hopf bifurcation point HB is located 
at cr i cr « 1.282. There are two CF bifurcations in Fig. 1, where stable limit cycles are 
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replaced by unstable ones. Between these two CF points, which occur at (Ji^cF 1 ~ 1-77 and 
&i,CF 2 ~ 1.83, two stable limit cycles coexist. Therefore, depending on initial conditions, one 
of the limit cycles will be selected in simulations of the glycolytic oscillator with substrate 
recycling. 

A general mechanism of turbulence in oscillatory reaction diffusion systems is wavenumber 
instability, i.e., instability of uniform oscillations against phase-like fluctuations [1]. In Eqn. 
(1-2), there are two different uniform oscillations that might undergo wavenumber instability. 
We want to provide evidence that these oscillations are stable against phase-like fluctuations 
for the parameters in Fig. 1. For the fast, uniform oscillations which originate from the Hopf 
bifurcation point shown by filled circles in Fig. 1, the stability condition can be obtained 
by reducing Eqn. (1-2) to the CGLE, 

A = (1 + ic )A - (1 + ic 2 )\A\ 2 A + (1 + zc x )AA (3) 

In Eqn. (3) A is the complex amplitude, and u, /3, 7 are real parameters. The CGLE has a 
uniform oscillatory solution, A = exp(i(co — Oz)t), that is stable if the condition 1 + C1C2 > 
holds. In the appendix, we calculated c , c\ and c 2 corresponding to Eqn. (1-2). Our results 
show that the uniform oscillations are stable for the parameters used in Fig. 1, for any 
D a > and £) 7 > 0. For D a = we find that C\ = 0. Hence, the parameter region 
we are interested in is deep inside the Benjamin-Feir stability region given by 1 + C\Ci > 0. 
Although the CGLE is valid only near the HB point, it is likely that the uniform oscillations 
will remain stable until the the next bifurcation in the system, i.e., CF 1 in Fig. 1 Q|. Next, 
consider the uniform oscillations with low frequencies. Unlike the case of fast oscillations, 
no analytic approach is available in this case. Note that oscillations shown by filled circles 
and diamonds in Fig. 1 occur at the same parameters. Therefore, it is rather unlikely that 
the slow oscillations undergo wavenumber instability, contrary to the fast ones. Thus, we 
can assume that uniform, slow oscillations are also stable. 

It is known that strong perturbations can switch oscillations from one stable orbit to an- 
other in the glycolyctic model with substrate recycling Q| . Therefore, even if both uniform 
oscillations in Eqn. (1-2) are stable against wavenumber instability, strong perturbations 
can excite the system by switching the oscillations. This kind of excitability, however, will 
not lead to turbulence; in the parameter interval [ct^cf 1 ! ^cf 2 ]; the fast oscillations will 
suppress the slow ones as time progresses. But for | cr^ — (t^cfA <C 1 where the fast oscil- 
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lations become unstable, it is apparent that a complete suppression of slow oscillations is 
impossible. Here, because of complex interactions between stable, slow and unstable, fast 
oscillations, interesting spatiotemporal dynamics might develop. Therefore, we carried out 
a detailed numerical study in the neighborhood of CF 1 . 

III. WEAK TURBULENCE IN A BIRHYTHMIC MEDIUM OF GLYCOLYTIC 
OSCILLATORS 

For numerical integrations of Eqn. (1-2) in one spatial dimension, we used the fourth 
order Runge-Kutta method. Diffusion terms were approximated by the finite difference 
method. Numerical parameters are Sx = 0.005cm, St = 0.05s. The system size is defined 
as I = NSx, where N is the number of spatial grid points. In this paper we present results 
for periodic boundary conditions, but we also tested the main results with no-flux boundary 
conditions. We also tested selected examples with smaller values of 5x and St for fixed /. 
Our simulations show that Eqn. (1-2) are sensitive to initial conditions. By choosing initial 
conditions as small perturbations of uniform, slow oscillations with large amplitudes, we 
found that these oscillations are stable for a { < a i>C F 2 - But, near and to the left of CF 1 , 
uniform, fast oscillations with small amplitudes are found to be unstable. They spiral out 
from unstable orbits towards the orbit of stable, large amplitude oscillations. For strong 
perturbations near the CF 1 bifurcation point, we found spatiotemporal irregular motions in 
Eqn. (1-2). 

Fig. 2 shows a gray scale plot of spatiotemporal dynamics in Eqn. (1-2). Oscillations 
between the white and black colors show large amplitude oscillations displayed by j(x,t). 
There are also oscillations with higher frequency and smaller amplitude in Fig. 2. Because 
the latter ones are unstable, they can not suppress large amplitude domains. Although 
uniform, large amplitude oscillations are stable against small fluctuations, phase slips created 
by strong initial perturbations cannot be eliminated as time progresses. As a result, spatially 
nonuniform distributions of concentrations are seen at given time moments, Fig. 3. On the 
phase plane, these nonuniform distributions generate motions attracted by unstable orbits 
around the inner cycle shown in Fig. 4. We found that such unstable orbits act as a weak, 
virtual heterogeneity emerging randomly. They cannot entrain the bulk oscillations, but in 
their presence, phase slips cannot be eliminated. Instead, persistent spatiotemporal irregular 
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motions develop. 

To characterize the irregular motions in Fig. 2, we calculated the maximum Lyapunov 
exponent X^p in 2N dimensional phase space First we made a very long run of Eqn. 
(1-2) to confirm that the turbulence is stationary. Then, by using the same initial conditions, 
we simulated Eqn. (1-2) and its linear system for computation of f° r 7\ = 2 • 10 5 s. We 
found that the largest Lyaponov exponent is positive and small, ~ 2 • 10~ 3 . We also 
calculated a two-point correlation function, C(x) =< -y(x ,t)'y(x + x,t) >, where < .. > 
stands for an average over space and time Q|. Fig. 5 shows that C(x) ~ const at small 
values of x, indicating strong local coupling and an absence of short waves. A power-law 
decay of the correlation function at intermediate values of x implies the presence of chaotic 
motions. We found that the slope is k « —0.15. We also found no significant variations 
of k and X^yap with changes of <Ji and I. The small values of n and X^™ > suggest 
that spatiotemporal irregular motions shown in Fig. 2-4 can be characterized as a weak 
turbulence. 

We found that in Eqn. (1-2), stationary irregular motions can develop only for certain 
initial conditions and system sizes. In simulations with different initial conditions and system 
sizes, we observed sudden collapses of turbulent dynamics. Collapse of turbulence in Eqn. 
(1-2) means a complete suppression of small amplitude oscillations. Thus, we defined the 
transient lifetime of turbulence t p as the time interval from initial conditions to the moment 
when all oscillators come within a distance d of the orbit of stable, slow, large amplitude 
oscillations. In our simulations we used d = 0.03. Following Ref. |l2j, we plot an average 
transient lifetime t p versus the system size I in Fig. 6. Here, each filled circles is an average 
of 20 simulations with different initial conditions. Fig. 6 shows that, as the system size 
increase, t p grows exponentially. 

For some initial conditions, when I is close to 2cm, the turbulent solution does not 
collapse. The inset in Fig. 6 shows the number of cases, among 20 different simulations, 
when a collapse of turbulence has not occurred by T = 10 6 s. (These cases were not included 
in calculations of the filled circles in Fig. 6.) We continued two cases in the inset (at 
/ = 1.75cm) up to T = 10 8 s and did not observe a collapse of motions near the inner cycle 
in Fig. 4. 

Numerical experiments indicate that if virtual heterogeneities reside sufficiently far from 
each other, a stationary pattern is possible in the interval o~i e [1.055, 1.075]. Fig. 7 gives 
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an example of such a pattern. Here virtual heterogeneities are located from each other by 
distances between 0.5cm and 1cm. Note that these quasi-periodic structures are not related 
to a Turing instability, which emerges due to differences in diffusion coefficients. Unstable 
oscillations at the onset of a CF bifurcation are the instabilities leading to these structures. 
The cellular structures in Fig. 7 are breathing because of the unstable oscillations. Numerical 
results show that as the system size increases, the cells breath coherently. 

A collapse of turbulence can be prevented if there is a reinjection mechanism for the 
unstable oscillations generated by the virtual heterogeneities. Naturally, wavenumber insta- 
bility can be such a mechanism. Using our calculations in the appendix, we simulated Eqn. 
(1-2) for parameters when the corresponding CGLE displays a wavenumber instability and 
found that small amplitude oscillations exhibit phase instability near Hopf bifurcation. For 
\a — <Ji ! cF 1 \ C 1 we found stationary weak turbulence. 



IV. WEAK TURBULENCE IN A CELL CYCLE MODEL 

In section III, we demonstrated that the CF 1 bifurcation point is crucial for turbulence 
in Eqn. (1-2). Mathematically, the term representing substrate recycling drives CF bi- 
furcations. In models of biochemical oscillations, terms representing enzymatic activities 
naturally arise. As an enzyme can quickly switch from being active to inactive and back 
again, ideal conditions for CF bifurcations exist in these models. Therefore, other biochemi- 
cal reaction diffusion models may also display the weak turbulence discussed in the previous 
section. As an example, consider a three variable model of the budding yeast cell cycle, 



dX 

— = m(h + k 2 T) - (k 3 + hY + k 5 Z)X + D X AX, (4) 

dY (fc 6 + fc 7 Z)(l-y) (k 8 m + k 9 X)Y a V f K\ 

-dt = J 1 + l-Y J 1 + Y + DyAY > (5) 

i rv 

— = (fc 10 + k n X) - k 12 Z + D Z AZ, (6) 

T = G(X,P,J 2 ,J 2 ), (7) 

G(a,b,c,d) — , (8) 

b — a + be + ad + J (b — a + bc + ad) 2 — Aad(b — a) 

where the transcription factor T for X is given by the Goldbeter-Koshland function G [13]. 
X, Y, Z are dimensionless variables and m is a dimensionless parameter. Our units of time 



and space are min and cm, respectfully. 

When Dx = Dy = Dz = 0, Eqn. (4-8) are the reduced version of a budding 
yeast cell cycle model [19]. Here, X represents the concentration of cyclin-dependent pro- 
tein kinase(CDK), Y and Z are concentrations of two different anaphase promoting com- 
plexes(APC), APC/Cdhl and APC/Cdc20 respectively. In Eqn. (4-5), m represents the 
cell's mass, which will be used as a primary bifurcation parameter. Eqn. (4-8) display CF 
bifurcations as shown in Fig. 8. For small m, Eqn. (4-8) also display saddle node bifurca- 
tions, a feature universal in cell cycle models |l4j,[l9j. Here, our concern is the neighborhood 
of CF 1 in Fig. 8. 

There are no experimental measurements of diffusion coefficients for CDK, and APC 
factors. But it is known that diffusion coefficients of average sized proteins in cytoplasm are 
order of lO" 4 ^ or smaller 2G, 211. As our goal is a demonstration of weak turbulence in a 

mm l " 1 "J ° 



representative model of biochemical oscillations, we choose Dx and Dz arbitrarily subject 
to this upper bond. For simplicity, we assume Dz = 0. 

For simulations of Eqn. (4-8) we used the same method as in the previous section with 
St = 0.05min, Sx = 0.005cm. We found numerically that for strong perturbations, Eqn. 
(4-8) display a weak turbulence, Fig. 9. Typically, for Dx < Dy, we found transient, 
weak turbulence. When Dx << Dy, numerical experiments lead to stationary turbulence. 
For instance, we simulated Eqn. (4-8) up to T = 10 7 min for m — 3, Dx = 6 • lO -7 ^-, 
Dy = 10~ 4 ^-, Dz = and I = 1.28cm and found stationary turbulence for a number of 
different initial conditions. 



V. DISCUSSION 

We have shown in this paper that two representative mathematical models of biochemical 
oscillations exhibiting birhythmicity, glycolytic and cell cycle models, display weak turbu- 
lence, intermiitency of large and small amplitude oscillations. We revealed that unsta- 
ble oscillations near cyclic fold bifurcations are the mechanism of transient turbulence in 
birhythmic media. In the presence of wavenumber instability, weak turbulence is stationary. 

Recently, Stich et. al. 22, l^] proposed an amplitude model for birhythmic media. 
An interesting question is whether the weak turbulence we discussed in this paper can be 
found in their model? First, let us mention two important differences between our models 



8 



and the amplitude model of birhythmic media. In our case, a cyclic fold bifurcation is 
crucial for turbulence, but the amplitude model describes a pitch-fork bifurcation of limit 
cycles. Secondly, both fast and slow oscillations in the amplitude equation are smooth, but 
in our case, slow oscillations are strongly anharmonic. Besides these differences, it is well 



known that if phase slips develop, the CGLE equation generates defects |24j . Thus, these 
facts indicate that instead of intermittency of small and large amplitude oscillations, defect 
turbulence is likely in the amplitude model of birhythmic media. 

To date, there are no experimental evidences of weak turbulence in glycolysis or in the 
cell cycle. Our results are pure theoretical predictions of mathematical models. The system 
sizes we simulated are much larger than the typical size of an yeast cell (10 _3 cm). Therefore, 
weak turbulence is not expected in yeasts. Interestingly, some slime molds grow as syncytial 
Plasmodia (many nuclei in a common cytoplasmic pool) that are many times larger than 



a typical yeast cell; cells 15cm in diameter can be grown in the laboratory |25j]. Waves of 
nuclear division are observed in these multinucleate plasmodia 2f], and, as we have 
shown it is possible that these waves exhibit weak turbulence. Note that weak turbulence 
in the cell cycle would mean irregular oscillations of CDK. But for a normal cell cycle, 
large amplitude oscillations of CDK are essential; CDK activity must drop below a certain 
threshold for nuclei to exit mitosis and divide. Therefore, hypothetically, weak turbulence 
in syncytial plasmodia might lead to mitotic arrest of certain nuclei in the plasmodium. 

A more quantitative characterization of wavenumber instability of unstable oscillations 
at the onset of a cyclic fold bifurcation, as well as simulations in two spatial dimensions [3] 
are problems in the future. 



APPENDIX A: COEFFICIENTS OF THE CGLE FOR A GLYCOLYTIC MODEL 
WITH SUBSTRATE INHIBITION 

In this appendix, following standard procedures in Ref. Q], we will calculate coefficients 
of CGLE for the glycolytic model. For a convenience we assume Q = 1 in Eqn. (1-3). First, 
let us find uniform steady state solutions ao and 70, 

70 = fx/ks, (Al) 

K\-2fi + a) + 7o 4 (-2(/i + a { ) + a) 

"o = z 

c 
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+ ALa~ab + (1 + 7o ) V 2 6 2 
(l + 7o)c 



(A2) 



where, a = if 4 /i+7Q(/i+<7j), 6 = (i^ 4 + 7 4 ), c = 2(K 4 (fj,— a)+ r y 4 (fi+(Ti — a)). Next we perform 
a linear stability analysis of (ao,7o) against small fluctuations 5a, <5 7 oc expiiqx + i\t). At 
the critical wavenumber q cr = 0, we obtain a characteristic equation, 



A 2 + (ai + a 2 + A: S )A + ai&; s = 0. 

In Eqn. (A3), a x and a 2 are given by 

a(l + 7o ) 2 (L + 2La + (1 + a ) 2 (l + 7o ) 2 ) 



(A3) 



ai 



a 2 



(L + (l + a )2(l+ 7o )2)2 
4K 4 7 ^ 2<TLa (l + a )(l+7o) 



(K 4 + 7o 4 ) 2 (L + (l + a ) 2 (l + 7 o) 2 ) 2 ' 
Let us define such a critical value for the bifurcation parameter cr, = <7j ;Cr . that 

Qi + a 2 + k s = 0. 



(A4) 
(A5) 

(A6) 



Eqn. (A6) is the condition for a Hopf bifurcation; the characteristic equation has pure 
imaginary solutions, Ao = ±.iy/a,ik s . 

Let a be defined by a — ?iz£i<£l_ We develop the Jacobian matrix L of Eqn. (1-3) in 

&i,cr 

powers of /i, 

L = L + fiLi + .... (A7) 



At \i = the Jacobian is given by 



-ai 



a 2 

-fc s - a 2 . 



(A8) 



We find the right u and left u * eigenvectors of L corresponding to A , 



u 



u ' 



1 



(A9) 



(A10) 
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We find further, 



Li 



4if 4 7gc^ cr 
(K* + 7o 4 ) 2 








-1 
1 



(All) 



Let us first find c in the CGLE. It is given by c = 7^x7, where 



Ai = u *Liu 



(K* + 7o 4 ) 2 

We see that Ai is pure real, therefore, Co = 0. Now following again [J, we find ci, 



(A12) 



D 











(A13) 



d! + zcf" = u *I>uo, 

Cl = d "/d! = [^(2^zJh 

Calculation of c 2 is a little more tedious. We need to find 



V+ = V_ = -(L - 2A ) 'MqUqUq, 
V = -2L " 1 MoU Uo, 



(A14) 
(A15) 

(A16) 
(A17) 



= 9 +W 



-2uo*M u Vo - 2u *M u V 
-3u o *N u o *uo*u . 



(A18) 



Parameter c 2 in the CGLE is given by a formula, c 2 = g" /g'. We find that c 2 = if , where 



</ = -3fc s [fc s m Q 2(2m tt 2 — m Q7 ) + ai((2m Q 2 — m ai ){m a 2 — m aa + m 7 2) 

— fc s (n Q 2 7 — 3n a 3)) + 3a 2 (n a7 2 — n a 2 7 + n a 3 — n^)], (A19) 

■[10/c 2 m Q 2 2 + aifc s (14m a 2 2 — 14m a 2m ai1 + m ai 2 + 10m Q 2m 7 2 + 9fc s n Q 3) + 

+a 2 (4(m Q 2 — m Qi7 + m^) 2 + 3& s (n a7 2 — 2n Q 2 7 + 3n a 3))]. (A20) 



In the above expressions, m c 



i i9 2 0(o,7) ' 



g a 2 ;«o,70' m «7 



' a 2 ^(q, 7 ) \ 

- 9a97 / Q 0,70' " l 7 
/ 9 3 0(q,7) \ 



(/^-4 +7 4-)3 I g 7 2 Mo, 70' ^a 3 V da 3 /Q!o,7o; ^a 2 7 " 

v da Q^ )a 0tl0 , n 7 3 = ( - g^3' 7 ^ )a ,7o- To save space we do not present here cumbersome 
expressions for these derivatives. 



11 



For parameters in Fig. 1 we find that <jj jCr = 1.282 and C2 ~ 2.21. From Eqn. (A15) we see 
that if D a = D y , c\ = 0. Therefore, I+C1C2 > for parameters used in this paper. If D a = 0, 
Eqn. (A15) gives the minimal value, c\ = —0.47. In this case 1 + C\C<i ~ —0.03, therefore, 
wavenumber instability is possible. However, turbulence must be weak as the parameters 
are very close to the stability condition 1 + C\C2 > jjj. A stronger wavenumber instability 
is possible, for example, for ll = 0.28s" 1 , K = 12, D a = 5 • 10~ 7 ^, £> 7 = 1 • lO" 5 ^ and 
other papameters are the same as in Fig. 1. In this case, we find that <7j jCr ~ 1.095 and 
1 + Cl c 2 = -0.416. 
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LIST OF FIGURES 

FIG. 1. A bifurcation diagram of Eqn. (1-3). HB marks a Hopf bifurcation point, CF 1 ' 2 
mark cyclic fold bifurcations. Parameters are: v = 0.25s -1 , n = 4, K = 11.5, a = lis -1 , 
k s = 0.05s -1 , L = 3400000. 

FIG. 2. Space-time pattern of 7 in a weak turbulent regime of Eqn. (1-2). The space and 
time spans are / = 1.75cm and T = 5 ■ 10 3 s. The pattern was obtained by recording 7(2;) 
with a time interval r = 5s. D a — D 1 — 1 • 10 -5£ ^- and a"; = 1.065s -1 . Other parameters 
are the same as in Fig. 1. 

FIG. 3. Snapshots of spatial distributions of a at two different time moments. Parameters 
are the same as in Fig. 2. 

FIG. 4. A phase plane view. The outer cycle shows the orbit of stable uniform oscillations 
with a period r = 300s. The inner cycle shows the orbit of small amplitude, fast oscillations 
with a period r = 290s at o~i = 1.08s -1 . With the decrease of <7;, the inner cycle disappears, 
but it still can attract neighboring trajectories creating a virtual, chaotic heterogeneity in 
Eqn. (1-2). The solid lines show oscillator distributions at two different time moments. 
Parameters are the same as in Fig. 2. 
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FIG. 5. A log-log plot of the spatial correlation function. Parameters are the same as in 
Fig. 2. 

FIG. 6 Average transient lifetime versus the system size. The inset shows the cases when 
a collapse of turbulence has not occurred by T = 10 6 s. Parameters are the same as in Fig. 
2. 

FIG. 7. Breathing periodic structures. / = 3.5cm, other parameters, as well as the time 
and space spans are the same as in Fig. 2. 

FIG. 8. Bifurcation diagram of a cell cycle model. Rate constants fcj are in units min -1 , 
h = 0.002, k 2 = 0.053, k 3 = 0.01, h = 2, k 5 = 0.05, k 6 = 0.04, k 7 = 1.5, k 8 = 0.19, 
k 9 = 0.64, k 10 = 0.005, k n = 0.07, k 12 = 0.08. Other parameters are P = 0.15, J x = 0.05, 
and J 2 = 0.01, I = 1.28cm, D x = 6 • 10~ 7 ^, D Y = 10~ 4 ^ and D z = 0. 

£. i mm ' 1 mm ^ 

FIG. 9. Turbulence in a cell cycle model. Space time plot of Y field in Eqn. (4-8). The 
space and time spans are L = 1.28cm and T = 2500min. The pattern was obtained by 
recording Y(x) with a time interval r = 5min. 
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